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Abstract 

Emulsion Stability Simulations (ESS) are used to estimate the coalescence time of one drop of hexadecane pressed by buoyancy 
against a planar water/hexadecane interface. In the present simulations the homophase is represented by a big drop of oil at 
least 500 times larger than the approaching drop (1/im to 10 fim). Both deformable and non-deformable drops are considered 
along with six different diffusion tensors. In each case van der Waals, electrostatic, steric and buoyancy forces are taken into 
account. The coalescence times are estimated as the average of 1000 random walks. It is found that the repulsive potential 
barrier has a significant influence in the results. The experimental data can only be reproduced assuming negligible repulsive 
barriers, as well as non-deformable drops that move with a combination of Stokes and Taylor tensors as they approach the 
interface. 
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I. INTRODUCTION 



Nowadays it is not yet possible to measure the coalescence time between two small drops of oil suspended in quiescent 
water. In an approximation to the problem four different experimental techniques are often used to determine: (a) 
the lifetime of thin films of macroscopic radii [1], (b) the nature of the forces between two fixed drops (Atomic Force 
Microscopy) [2, 3], (c) the dynamic interaction forces between two particles during flow [4], and (d) the coalescence 
time of oil droplets pressed against a planar oil- water interface [5-7]. The first technique allowed identifying the 
stages of the coalescence process between deformable interfaces. The second showed significant differences between 
the behavior of films and droplets during coalescence. The third method outlined the importance of the Brownian 
movement in the scattering of particles of micron size covered by /3-casein. The fourth technique evidenced an 
interesting dependence of the coalescence time with the radius of the approaching drop. 

Using optical microscopy Dickinson et al. [6] reported the coalescence time of micron-size (1 to 5 /zm) droplets of 
hexadecane pressed by buoyancy against a planar interface. For this purpose, a small drop of oil was released in an 
aqueous protein solution, just below a planar water/hexadecane interface. The coalescence time between the droplet 
and the interface was determined by direct observation using optical microscopy. Proteins (/3-casein, /-c-casein and 
lysozyme) were used in order to guarantee the immobility of the surfactant layer at the oil/water interface. However, 
the surfactant concentration was kept low enough in order to guarantee the coalescence of the drops. It was found 
that independently of the protein used, the coalescence time decreased with the increase of the radius of the droplets. 

Ivanov and Kralchevsky [8] studied the hydrodynamic flow between a sphere and a planar interface, emphasizing on 
the effect of the drop deformability on the coalescence time. They predicted a minimum of the coalescence time as a 
function of the drop radius for the case of tangentially immobile interfaces and a fixed driving force. Basheva et al. [7] 
studied the asymptotic limits of this size dependence both theoretically and experimentally. Using an experimental 
set up similar to the one of Dickinson et al., they reported a decrease of the coalescence time of soybean oil drops 
(2 ^m < r-j < 100 /im) pressed against a planar water/oil interface. For these experiments, Bovine Serum Albumin 
(BSA) was used as a surfactant. 

In order to explain their experimental results, Basheva et al. [7] argued that when a drop is released in a bulk liquid, 
it moves according to Stokes law until it approaches the interface to a sufficiently small distance. Then it markedly 
decelerates due to the increase of the viscous friction in the remaining gap. Hence, the authors defined the lifetime of 
the drop at the interface, r, as the time elapsed between the moment in which the drop start to move slowly until it 
disappears by coalescing with the large homophase: 



(1) 



t v{hy 

where v(h) is the velocity of thinning of the liquid film between the drop and the interface, h is the closest distance of 
approach between their surfaces, h m i is the gap width at which the thinning begins, and /i cr i t , the minimum distance 
that can be attained (critical thickness) before the film breaks and coalescence occurs. If the interfaces of the drop 
and the planar boundary are tangentially immobile, and if the drop keeps its spherical shape until contact, then the 
expression of Taylor for the drainage of the intervening film between hard spheres can be used [9] 
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Here F is the external driving force, r\ is the dynamic viscosity of the external liquid, is the radius of the small 
droplet and rj is the radios of the large drop. Making rj — » oo: 



hF 



Basheva et al. assumed that the velocity of the approaching drop could be expressed as a combination of the 
velocities predicted by the tensors of Taylor and Stokes: 

- = — + , (4) 
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where vm is given by Eq. (3) and vst — F/6irr]ri refers to the Stokes law for motion of a sphere in an unbounded 
liquid. 
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Substituting Eq. (4) in Eq. (1), and using the expression for the buoyancy force: 



F = ^irrfApg, (5) 

where Ap is the density difference and g is the gravity, an inverse dependence between the lifetime of the drop at the 
planar interface, and its radius is obtained: 
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According to Eq. (6), the lifetime diminishes with the increase of the droplet radius r^. This trend is just the 
opposite to the one found for drops of hundreds of microns (millimeter size) . In the latter case the drops deform close 
to the interface, producing a circular film. This film drains according to Reynolds law, producing an increase in the 
lifetime of the drops as a function of their radii [7] . 

It is remarkable that the experimental data on the coalescence of a drop with a planar interface follows Taylor law, 
since this expression was initially deduced for the case of two colliding spheres. However, it is well known that the 
radius of curvature of a drop increases with its size, approaching asymptotically the radius of curvature of a planar 
interface. Conversely, the behavior of the drop/interface system could be simulated using two drops of very distinct 
sizes. One big drop fixed in space and sufficiently large to be regarded as a planar interface by the approaching 
droplet, and one moving drop of micrometer size rising as a consequence of the buoyancy force. 

In this work the algorithm of Emulsion Stability Simulations (ESS) is used to reproduce the data of Dickinson 
et at regarding drops between 1 and 10 /im. It will be shown that the simulations fit the experimental data if the 
approaching drop is assumed to be a non-deformable sphere, moving with the diffusion tensors of Stokes or Taylor 
depending on the distance of approach. It is also demonstrated that the use of a truncated spheres as a model of 
dcformable drops does not fit the experimental data for these range of sizes, independently of the tensor employed. 

The article is structured as follows: In Sec. II an overview of Emulsion Stability Simulations is presented, in Sec. 
Ill a description of the calculations is given along with some technical details, in Sec. IV the paramctrization of the 
potentials is explained, and some preliminary results concerning the form of the potentials are presented. In Sec. V 
we show the results of the calculations and Sec. VI presents the conclusions. 



II. EMULSION STABILITY SIMULATIONS 



Emulsion Stability Simulations are based on the algorithm of Brownian Dynamics published by Ermak and Mc- 
Cammon [10-12]. If the divergence of the diffusion tensor is assumed to be negligible, the displacement of particle i, 
Yi(t + At) — Yi(t) is equal to: 



r l {t + At) = r t (t) + ^^At + H, (7) 

where the second term on the right hand side accounts for the effect of conservative forces on the particle movement, 
Di is an effective diffusion constant (tensor) of particle i, is the total force acting on i, k is the Boltzmann constant, 
T the temperature, At the time step, and R is a random term representing the Brownian motion of the particle. 

The diffusion constant is equal to Di = D fcolvfcarr where D a — kT/Q-K-qri is the diffusion constant of Stokes. The 
first correction term /iorr, takes into account those factors that change the expression of the diffusion constant at 

(2) 

infinite dilution [13]. The second correction term / C or r takes into account the hydrodynamic interactions between the 
particles caused by the movement of the surrounding liquid as the particles advance [11]. 

In the most common case, a calculation begins distributing a set of oil drops in a cubic box of side length L. It 
is assumed that the molecules of oil mainly determine the van der Waals interaction between the particles. Instead, 
the repulsive interactions depend on the amount and chemical nature of the surfactant molecules adsorbed to the 
interface of the drops. The program has several routines for apportioning surfactant molecules amongst the drops. 
Once the surfactant has been distributed, the surface properties of the drops (like charge, intcrfacial tension, etc.) 
can be calculated. Then, the diffusion constant and interaction forces can be computed and the drops moved using 
Eq. (7). At every time step, the program checks for the coalescence of drops. In the case of non-deformable drops, 
coalescence occurs whenever the distance of separation between the centers of mass, r^ , is smaller than the sum of 
the radii of the drops. When this happens, a new drop is created at the center of mass of the coalescing particles. 
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The present version of the code can either simulate the behavior of non-dcformablc or dcformablc drops. In both 
cases the particles follow the same equation of motion, Eq. (7), but the analytical form of the diffusion tensor and 
the interaction forces change. 

If the mode of deformable droplets is selected, it is assumed that the deformation of the drops occurs independently 
of the energy required for this process. Due to its simplicity, the model of truncated spheres is used [14-16]. According 
to this model, three regions of approach between two particles are defined: 

Region I: If the distance of separation between the centers of mass of the drops, , is larger than rt + rj + h[ n i, where 
/lini stands for the initial distance of deformation, the drops maintain their spherical shape. Consequently, the 
radius of the liquid film between the flocculating drops is zero, rf = 0. 

Region II: The drops change their shape from spheres to truncated spheroids. This region covers the range of 
distances between the beginning of the deformation rf ^ 0, and the attainment of the maximum film radius: 
r i = r fmax = V r ihini- In this zone, the close distance of separation between the surfaces of the drops is assumed 
to be constant h = h m i [14], and: 
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(8) 
(9) 



Region III: The film radius already attained its maximum value, n — n max = V^ini, and the intervening liquid 
drains until it takes a critical distance of approach: 
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Accurate estimation of the initial distance of deformation hi n i, is very difficult since it results from a balance between 
hydrodynamic and interaction forces. The combination of the movement of the particles with the numerical solution 
of the exact formula of ft.; n i [16] is too demanding in terms of computational resources. Hence, we fitted the curves 
obtained in Rcf. [16] with a polynomial expression for the approximate estimation of h- ln i(ri, 7): 



hi„i = [1.2932 x 10 8 - 8.6475 x 1CT 9 
x exp(-ri/1.8222 x 1CT 6 ] 

3.3253 + 5.9804 exp(- 7 /0.00402) 
X 3.3253 + 5.9804 exp(- 10" 3 /0.00402) ' ( ' 

where 7 is the interfacial tension. For the value of fr C rit, the expression published by Scheludko and others is used 
[17-20] 



A A \V 4 



»-n-5FM • (13) 

where A cr ; t — n/10 and Ah is the Hamaker constant. 

Even in the mode of deformable drops, the particles behave as spheres if > + rj + /i; n ;. This means that 
the potential of interaction and diffusion constant correspond to the ones of spherical particles within Region I. At 
h = hi n i, the code calculates the dimensions of truncated spheres which are compatible with the actual distance of 
separation between the centers of mass of the spherical drops (r^ < + rj + /ij n i)- In this case, the expressions of the 
potentials corresponding to two truncated spheres are employed. Different expressions for the interaction potentials 
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of truncated spheres are available from the literature [16]. They include, van der Waals, electrostatic, etc, see Table I. 
Those potentials are expressed in terms of the width of the film and the particle radius. Use of Eqs. (8)-(ll), allows 
algebraic differentiation of the potentials in terms of . As a result it is possible to obtain analytical expressions for 
the force. Those expressions can be evaluated using the value of the film width and film radius corresponding to each 
region of approach. 

It is important to remark that a pair of drops does not necessarily move sequentially between regions I, II and III. 
The movement of the particles described by Eq. (7) may lead a couple of particles to go into the deformation region 
and back as a consequence of the random movement of the particles and/or their interaction forces. The program uses 
Eq. (7) to move each particle separately, although the analytical form of the forces and diffusion constants depend 
on the relative distance between the particles. 

Table I shows the analytical form of the potentials employed in the present calculations. Notice that the geometrical 
part of the van der Waals and electrostatic potential changes as a function of deformation. Additionally, two new 
potentials appear during the evolution of the film (Region II). They take into account: (a) the surface deformation 
energy coming from the increase of interfacial area as the spherical drops turn into truncated spheres, and (b) the 
bending elasticity potential related to the curvature of the interface [14]. These two potentials change with the 
interparticle distance during the formation of the film (Region II) , but reaches a constant value once a maximum film 
radius has been attained. Hence, they do not contribute to the value of the force in Region III where Eqs. (10-11) 
hold. 

Notice that in the present calculations the same form of the steric potential for spherical and deformable droplets 
is used (Table I). This was done on purpose for several reasons. First, it avoids some anomalies in the behavior of the 
potential that might occur when the volume of overlap changes as a consequence of the sudden transition between 
spherical particles and truncated spheroids (Regions I and II). Second, it simplifies the parametrization of the steric 
potential produced by casein when it is adsorbed to a liquid interface [4, 21-23]. Third, structurally different proteins 
like lysozyme, ft-casein and /3-casein showed a similar behavior of the coalescence time as a function of the particle 
radius (for a planar interface aged during 20 minutes) [6] . Therefore the exact analytical form of the steric potential 
does not appear to be very significant. Fourth, according to the experimental methodology, bare drops of hexadecane 
were formed by vigorous mixing with a buffer solution. Then the drops were injected into the cell containing the 
protein solution and the planar interface. All the drops reached the interface within the first five minutes, and 
frequently before 30 s. However, it takes hours for /3-casein to attain an equilibrium concentration in the presence of 
a hexadecane/ water interface (see Fig. 9 in Ref. [6]). Consequently, the amount of protein adsorbed at the interface 
of the drops is uncertain, and the exact magnitude of the steric barrier is unknown. As a result of all these factors, 
it was considered convenient to simplify the form of the potential as much as possible, and test the effect of different 
steric barriers on the outcome of the simulations. 

Fig. 1 illustrates the methodology employed for evaluating the diffusion tensor of the drops. In the case of spherical 
drops the space around each drop i is divided into two regions. An imaginary sphere of radius di n t = 2r^ delimits the 
internal region. If a neighbor particle reaches the internal region of particle i (case I): h < ri (with h = rij — ri — rj), 
the position of the closest particle is used to calculate the diffusion constant of i. The program has the possibility 
to select between the expressions of Honig [29] , Taylor [9] , and a linear combination of the expressions of Taylor for 
immobile and mobile interfaces. If none of the surrounding particles reaches the internal region of particle i (case 
II), the expressions of Stokes, Mills and Snabre [30], Honig et al. [29], Taylor [9], or Beenakker et al. [31-33], can 
be assigned to particle i. The formulas of Mills et al. and Beenakker et al. use the volume fraction of oil in the 
simulation box to evaluate an empirical form of the diffusion tensor. 




(a) (b) 



FIG. 1: Calculation of the diffusion tensor, (a) Spherical: Inner region (0 < h < n), outermost region (h > n), h = r»j —ri — rj. 
(b) Deformable: Inner region (0 < h < /imi), intermediate region (hi n i < h < n), outermost region (h > n). 

In the case of deformable drops the options of Stokes, Mills et al, Beenakker et al, and Honig et al. are also 
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TABLE I: Potentials used in the simulations. In these equations, n is the radius of the small droplet, rj is the radius of the 
large drop and ft is the minimum distance between their surfaces. For the van der Waals (vdW) potential: Ah is the Hamaker 

/ = ft + n + \J r\ — r 2 , L = rj + . r 2 — r 2 , d = \JK 2 + 4r 2 , ft and rt are the thickness and radius 



constant, x 



h 
2r 4 ' 



V = 



-Coi, z is the charge number, e is the dielectric 



of the film, respectively. For the electrostatic potential (elect): k 2 = - 
permittivity of the medium, C e \ is the electrolyte concentration, e is the electron charge, kT is the thermal energy, *& s i and * s j 
are the surface potentials for the small and large drops, respectively. For the steric potential (st): V v is the molar volume of 
the solvent, \ is the Flory-Huggins solvency parameter of the protein, 4>j and <f>i are the average volume fraction of the protein 
around each sphere, <pi = — ^ ™ X \3 — gr, with r the number of molecules per unit area, p p the density of the protein, M p 

the molecular weight of the protein and 8 the width of the protein layer. Volumes v a , Vb, and v c depend on h. Their explicit 
geometrical expressions to calculate the volume of overlap between the interacting particles can be seen in Ref. [26]. For 



the dilatational (extensional) potential (dil): 70 is the interfacial tension and r a 
B = 1.6 x 1CT 12 N [27, 28] is the interfacial bending moment. 



2r i T 3 
ri+rj 



For the bending potential (bend): 



available if all the neighbor particles are located within the external region of drop i. However, the internal region 
around each particle is now subdivided into two parts. The outermost zone of the internal region corresponds to 
Region I defined above in terms of the approximation distance between two deformable drops. In this region the 
deformable drops are close to one another but still maintain their spherical shape, h[ D i < h < ri. Consequently the 
expression of Honig et al. is used to evaluate the diffusion tensor. The innermost zone of the internal region of particle 
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i corresponds to the regions of approach II and III defined above. Within these regions, the drops behave as truncated 
spheres. Hence, the set of formulas compiled by Gurkov and Basheva for deformable drops (see Table 1 in Ref. [34]), 
can be selected. 

It might occur that the range of validity of the mathematical expression used to simulate the diffusion tensor at 
close separation distances is exceeded in some cases. The program checks if either fcolr > 1 or /iorr > 1. In that case, 
it has the option to correct the diffusion tensor using the expression of Honig et al, or use the calculated diffusion 
constant regardless of its value. 

The analytical form of the set of tensors used in the present simulations are shown in Table II (see the section on 
Computational Details). 
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Stokes Immobile: A sphere in an 
unbounded liquid [7]. 


Dqt = Dn 


II 


Honig [35]. 


Dh " = 6^+13^+2 D °- 


III 


Taylor Immobile: Two spheres of radii 
ri,rj, and h <C n, rj [34]. 


D T , = 4D ^h. 


IV 


Taylor Mobile: Two sphererical liquid 
drops or radii ri, rj [34]. 


7-1 ^j-i ^ u ( i-i.7ne-o.46ie 2 A • 
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Reynolds Immobile: Film between two 
circular disks or radius n and h -C r{ [34]. 


DRe=4nD !£. 


VI 


Two deformed drops [15]. 





TABLE II: Tensors used in the simulation. Here n is the radius of the small droplet, rj is the radius of the large drop, 
r* = h and n are the thickness and radius of the film, Do = 6 ^ r . , £ = ^"\/?' ^ s ^he dynamic viscosity of the 

continuous phase, 77, is the dynamic viscosity of the disperse phase, u — r,J ^ Tj , where rij — n + rj + d is the distance 
between centers with Rq as a radius of reference and es have values between 0. 0&1 and 1. 

While the coalescence of non-dcformable drops occurs when the particles overlap, deformable drops could coalesce 
through two different mechanisms: 

1. The drainage of the intervening film between flocculated drops. In this case the drops move until the distance 
of separation between their surfaces reaches h — h cr i t - 

2. Surface oscillations and/or the formation of holes promote the rupture of the film prior to its drainage [1, 17, 19, 
36-40]. In this case, the surface oscillations of the film are not simulated explicitly, but a stochastic procedure 
is implemented, taking into account some relevant aspects of the rupture process. 

First, the lifetime of a doublet, (t^), is calculated using one out of two procedures: 

2.1. Counting the time continuously from the moment a doublet enters regions II or III until either h — Z^crit 
(coalescence occurs) or the doublet separates: > + rj + /i ini . 

2.2. Adding the time steps in which each couple of particles enter the regions II or III. In this case, r%j is 
different from zero after a pair of particles enter regions II or III the first time. This way of counting the 
time assumes that the probability of rupture is proportional to the time spent by the doublet in the region 
of deformation. 

Following, a random number between —1.0 and 1.0 is assigned to the surface of each drop. The amplitude of 
each capillary wave (Ai) is estimated as the product of the referred random number times the value of h cr i t - 



A = ran(t) x h CTit . (14) 
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If the mechanism of surface oscillations is activated, the value of 77,- is compared at each time step with a 
characteristic time deduced by Vrij [38] for the fastest increase of surface oscillations: 



rvrij = 967^777/4; A H 2 . (15) 

Eq. (15) was deduced assuming van der Waals interactions only. A more general expression for TVrij [38] 
requires knowledge of the second order differential of the free energy in terms of h, under certain mathematical 
restrictions. Hence, it is very difficult to calculate when the model of dcformable drops is employed. Therefore 
the program has the option to introduce the value of ryrij as part of the input data. 

The value of the tension in Eq. (15) is approximated by the average of the interfacial tension of the two drops 
(7»+7j)/2. The value of ji is then calculated at each time step of the simulation, using the number of surfactant 
molecules adsorbed: 



7* = 7o + (7cmc - 70 ) yj^tk J ■ ( 16 ) 

Here 70 and 7 cmc stand for the value of the O/W interfacial tension in the absence of surfactant molecules, 
and at the CMC of the surfactant employed. -/V Sj , being the number of surfactant molecules adsorbed to the 
interface of drop i, and N™ ax , the maximum number of surfactants that can be adsorbed to that drop. When the 
surfactant concentration is equal to zero, N s ^=0 at all times, and 7* = 70. In the case in which the surfactant 
concentration is enough to cover the drops completely, the equilibrium surface tension 7$ = 7 cmc is reached. 

Coalescence occurs whenever the total height of the surface oscillations is greater than h cr it- The height of the 
oscillations is approximated by: 



exp 



\ TVrij J 



(17) 



Eq. (17) takes into account that: (a) the surface oscillations increase exponentially with time [36-38], and (b) 
the capillary waves can be in-phase or out-of-phase. 



III. COMPUTATIONAL DETAILS 



In order to reproduce the lifetime of a drop of hexadecane at the water/hexadecane interface, two drops of very 
different sizes are used. The radius of the small droplet r (r"j = r), was varied between 1 to 10 /im. A fixed drop 
of 500 /im represents the interface. A square simulation box with a side length L equal to 4 R + d (where d is the 
initial distance between the drops, and R (rj — R) is the radius of the large drop) was employed. Figure 2 illustrates 
the proportion between the big drop and the moving droplet, as well as the spatial position of the doublet at the 
beginning simulation. The mean coalescence time and its standard deviation was estimated using 1000 random walks 
for each particle size. For each individual simulation the initial values of the random term in Eq. (7), R, were changed 
in order to favor different trajectories. The logistic equation in the region of chaos was used to avoid the saturation 
of the random number generator (see Ref. [35] for details). 

The surfactant concentration of /3-casein and the ionic strength of the aqueous solution employed in the simulations 
correspond to the ones used in the experimental measurements of Dickinson et al. (10~ 4 wt% protein, pH=7.0, ionic 
strength 0.1 M, T=298K [6]). 

As pointed out by Damodaran and Rao [41], the interfacial areas of the proteins resulting from T vs. t or n vs. t data 
are considerably smaller than the cross sectional areas of the native or unfolded proteins, which typically fall in the 
range 100 — 1000 nm 2 [41]. According to the experimental data of Graham and Phillips [42-44], the maximum surface 
coverage of /3-casein at the air/water and the oil/water interface occurs between 2 — 3mg/m 2 . Chen and Dickinson, 
determined the value of T for a hexadecane-in-water emulsion stabilized with 0.4% wt /3-casein, obtaining a value of 
2.95 mg/m 2 [45, 46]. This value suggests an interfacial area of the order of 14 nm 2 . Conversely, one can estimate the 
interfacial area of a /3-casein micelle from its radius (13 nm [47, 48]), and compute the number of molecules in the 
micelle from its molecular weight (250 kDa [47, 49]). In this case, a much larger interfacial area of 204 nm 2 is obtained. 
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Oil droplet 
of radius r 



2R 



FIG. 2: Model of the drop/interface system employed in the simulations. Here, r is the radius of the small droplet, and R the 
radius of the large drop resembling the interface, d is the initial distance between the surfaces of the drops. 



This discrepancy can be partially due to the internal structure of the casein micelles [50, 51], and the large number 
of conformations available for a long polyelectrolyte when it is partially dissolved in a liquid medium. 

The analytical form of tensors employed in the present simulations, D corr , are shown in the Table II. Following 
the scheme of Fig. 1, two tensors must be defined for spherical and deformable droplets, respectively. The first one 
corresponds to the outer region of approach denoted by r oxt in Fig. 1. The second one corresponds to the inner region 
of approach (r; nt in Fig. 1). 

In the case of spherical droplets the following combinations of tensors Axt/Ant, were tested: 

A. Stokes/Honig. Tensors 7 and 77 of Table II are used. The radius of reference i?o for the equation of Honig et al. 

was approximated by the average radius of the drops: 7?o = {fi + rj)/2. 

B. Honig/Honig. Tensors 7 and 77 of Table II are used with R = (r^ + rj)/2. 

C. Honig/Honig. Tensor 77 of Table II is used in both the internal and external regions of drop i, with Rq = r», 

where r t is the radius of the smallest drop. 

D. Stokes/Taylor Immobile. Tensors I and III of Table II are used. 

E. Stokes- Taylor Immobile/Taylor Immobile. Tensors I and III of Table II are used. The tensor of Taylor is used 

(2) 

in the external region as long as / CO rr = D^/Dq < 1. Otherwise, the expression of Stokes is used. 

(2) 

F. Taylor Immobile/Taylor Immobile. Tensor III of Table II is used regardless of the value of / CO rr- 

G. Stokes/Taylor Mobile. Tensors I and IV of Table II are employed. Values of 3.032 x 10~ 3 Pas and 8.905 x 10~ 4 Pas 

were used for the viscosity of hexadecane [52, 53] and water at T = 298.15 K, respectively. 

H. Stokes/Taylor Mobile- Taylor Immobile. Tensor I of Table II is used for the external region. Within the internal 

region a linear combination of tensors 777 and IV are used: 



A = C TR D Ta + (1.0 - C T R)D Ta , m , (18) 

where Ctr = N s s/N™ ax . It is clear from Eq. (18), that the total coverage of a drop generates the use 
of the Taylor's formula for immobile interfaces. However, it is known that the adsorption of /3-casein to an 
hexadecane/ water interface is very slow [6]. Hence, it is likely that in the experiments of Dickinson et al, the 
protein layer at the surface of the emerging drop is very dilute even at the moment of the collision with the 
planar interface. This favors an intermediate mobility of liquid at its interface, similar to the one described by 
Eq. (18). In the present simulations the temporal dependence of the protein adsorption was not considered. 
Hence, an intermediate interfacial mobility Eq. (18) can only be produced if either the amount of protein 
molecules dissolved in the aqueous phase is subtracted from the total surfactant population, or if a lower protein 
concentration is arbitrarily introduced in the simulations. Consequently, in order to study the effect of tensor 
H on the coalescence time of hexadecane droplets, a protein concentration of C s = 10 -9 M was used. This 
concentration generated a value of Ctr, = 0.33. 
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In the case of deformable droplets, the inner and outer regions of approach are separated by an intermediate zone. 
In this zone the approaching drops are still spherical. As shown in Fig. 1, the program always uses the expression 
of Honig et al. in this intermediate region. Notice also that the expression of Reynolds (Eq. V in Table II) requires 
h <C ri,rj. However, there might be cases in which the actual situation falls out of the range of validity of Eq. V. 
This might occur if the initial deformation distance is introduced as part of the input, or when the motion of particles 
of nanometer size is simulated. In order to account for this deficiency, it is assumed that the formation of a plane 
parallel film will necessarily delay the coalescence process. Hence, the code uses Eq. V as long as Z}R e /A) < 1. 
Otherwise, the expression of Reynolds is substituted by the one of Honig et al. in the internal region. 

The following combinations of Dext/Mnt were tested for the case of deformable droplets: 

I. Stokes/Reynolds. Tensors I and V of Table II are used. The value of hi n i is calculated by the program using Eq. 
(11). 

J. Honig/Reynolds. Tensors II and V are used. The value of /i; n i is calculated by the program using Eq. (11). 

K. Honig/Reynolds. Tensors II and V are used. The value of h- ln i is set equal to 20 /im and the value of rf = r,/10. 

L. Honig/Reynolds. Tensor II and V are used. The value of /i; n ; is set equal to 20 /im and the value of rf = r/100. 

M. Stokes/Danov. Tensors / and VI are employed. The value of h m { was calculated by the program using Eq. (11). 
Values for es = 1 and 0.1 were tested. 

For every combination of tensors a recursive calculation of the coalescence time was executed 1000 times. The results 
were compared with the experimental coalescence times obtained by Dickinson et al. These times were obtained by 
digitalization of Fig. 4 of Ref. [6] using the Engauge Digitizer. This data is represented by stars in Figs. 14-17. 

Most of the parameters employed in the simulations are shown in Table III. The mass of a casein protein was 
taken from Mobius: 24kDa [54]. The interfacial area was calculated using the radius (13 nm) and weight (250 kDa) 
of a spherical /3-casein micelle [47, 48] . The effective charge of the protein was evaluated varying the surface charge 
of a small hexadecane drop covered with /3-casein, in order to reproduce its surface potential £=— 27.8mV (<r = 
0.002161 C/m 2 ) [47]. Notice that the effective charge of the protein shown in Table III: 2.864e-, is substantially lower 
than the value reported by Dickinson 15e- [55] for the same molecule at pH= 7.0. In previous works, the described 
procedure generated an effective charge of 0.21e- for every single formal valence of a surfactant molecule. The value 
obtained for casein is 13.8 times higher than a single effective charge. This is reasonable in the light of our previous 
results [11, 56]. 



Hamaker constant [57] 4.90 x 1CT 21 J 

Ionic strength [6] 0.1 mol/1 

Surf, concentration [6] 4.16x 10~ 8 mol/1 
Molecular mass [54] 24 kDa 

Specific volume [58] 0.743 cm 3 /g 
Electric charge — 2.864e 



TABLE III: Parameters of the one /3-casein molecule. 

It is believed that the stability of milk in the presence of large salt concentration is related to the steric barrier 
produced by its casein proteins. However, that surface activity is mostly ascribed to K-cascin. This molecule contains 
63 hydrophilic aminoacids which lie on the outside of the casein micelle [50] . In the case of /3-casein only 50 out of 209 
amino acids are hydrophilic whereas the rest are mainly hydrophobic [48]. Depending on the protein concentration 
and the temperature, the core density of the /3-casein micelles varies between 0.4 and 0.9g/cm 3 , whereas the density 
of a globular protein is close to 1.35 g/cm 3 . Moreover, the density of the outer shell of the micelles is reported to be 
much lower, between 0.025 g/cm 3 and 0.14g/cm 3 [48]. 

In general, the steric potential of a protein is difficult to parametrize. V^t (Table I) requires the width of the steric 
layer (<$), the value of the Flory-Huggins interaction parameter (x), and the volume fraction of protein in the steric 
layer. 

Dickinson et al. [59] studied the adsorption of /3-casein at air/water and oil/water interfaces using Neutron Reflec- 
tivity. According to those measurements, the distribution of protein normal to the interface is well described by a 
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dense inner layer of 2 nm thickness located directly at the interface, and a more tenuous secondary layer of thickness 
5-7 nm extending into the aqueous phase. Hence, we selected a conservative value of <5=6.2nm for the width of the 
polymer layer outside the drop. 

The theory of Scheutjens-Fleer was used in a previous theoretical study regarding adsorbed /3-casein [21]. That 
formalism ascribes different values of the Flory-Huggins parameters to different segments of the molecule, depending 
on their chemical nature. We did not find either an experimental determination or a theoretical evaluation of \ for 
the whole /3-casein. It can be observed in Table I, that the difference (1/2 — x) is directly proportional to the value 
of the steric potential. Hence, x nas to t> e lower than 0.5 in order produce a repulsive barrier (0 < x < 0-5)- Tuinier 
and Kruif used the theoretical expression of Alexander-deGennes [50] in order to account for the steric interaction 
between K-casein micelles. According to that research, the height of the brush of protein is hardly affected by the 
solvent quality above pH 4.2. However, the relative height of the layer with respect to its value at pH 7, changes 
between and 0.2 at pH 2 when x is equal to 0.5 and 0, respectively. A value of x = 0.4 produces an intermediate 
relative height of 0.1 for the brush at low pHs. In the absence of other guideline, we selected a value of 0.4 for the 
Flory-Huggins parameter. 

In the ESS program the value of the volume fraction of protein in the steric layer around a drop can either be 
introduced as an input of the simulation or calculated. As shown in Table I, there is a simple equation that relates 
the volume fraction of protein around a drop (pi, with the number of molecules per unit area adsorbed to its interface: 



$. ffl^ (19) 

p P N A [(r i+ Sr-rfY liyj 

Use of the specific volume reported by McMceking [58]: 0.743 cm 3 /g for a casein molecule in solution along with 
an interfacial area of 204 nm 2 , produces the solid line illustrated in Fig. 3 for the interaction between two K-casein 
micelles (4> p = 2.1 x 10~ 2 ). 
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FIG. 3: Steric potential between two 100 nm k— casein micelles. Solid line: V B t (4> P = 2.1 x 10~ ); dashed line: V e t (<j> P = 
1.0 x 10~ 2 ); dotted line: Steric potential reported by Tuinier and Kruif (Fig. 2 in Ref. [50, 51]); dot-dashed line: V s t 
(4> p = 4.7 x 10 -3 ); dashed double-dots line: Vet (4> P = 7.9 x 10~ 4 ). 

Figure 4 illustrates the total interaction potential between a drop of hexadecane of 10 ^.m and the 500 fim drop 
representing the planar interface. The values of V^t were calculated using an interfacial area per protein molecule of 
204 nm 2 , and a volume fraction of protein equal to (f> p = 2.1 x 10~ 2 . The resulting repulsive barriers are enormous. 
These large barriers will completely prevent the coalescence of the drop of oil with the planar interface, in contradiction 
with the experimental evidence. Much larger steric barriers are obtained if the interfacial area of the protein is 
approximated by 14 nm 2 . 

It is known from the experiments of Dickinson et al. [6] that all emulsion drops coalesce with the planar interface 
as long as the interface is aged for only 20 minutes. In previous papers our group demonstrated that barrier heights 
AV larger than 12.7 kT prevent the coalescence of spherical drops [26, 35]. The barrier heights were estimated as 
the difference between the values of the maximum and the secondary minimum of the total interaction potential 

(AV = V max - V m i n ). 

It is possible to obtain surmountable repulsive barriers for the total potential varying several parameters of the 
steric potential. The magnitude of the electrostatic contribution is lower than the steric interaction due to the high 
value of the ionic strength (I = 0.1 M). However, the parameters of the steric potential are not independent. For 
instance, the volume fraction of protein around the drops depends on the width of the protein layer, which in turn 
depends on the surface excess of the protein at the interface (r). Due to the experimental set up [6], the value of 



11 



14000 
12000 
1 0000 

I— 

S 8000 
o 

■| 6000 
~ 4000 
2000 


-2000 ■ 



6 8 
d (nm) 



FIG. 4: Total potential of interaction between a large droplet of 500 /im and a small droplet of micron-size. For this calculation 
the volume fraction of protein around the drops was assumed to be equal to ^ p = 2.1 x 10~ 2 . Solid line: r = 2 pm; dashed 
line: r = 4 pm; dotted line: r = 6 pm; dot-dashed line: r = 8 pm; dashed double-dots line: r — 10 /im. 
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FIG. 5: Total potential of interaction between a large droplet of 500 /im and a small droplet of micron-size. For this calculation 
the volume fraction of protein around the drops was assumed to be equal to (j> p = 1.0 x 10~ 2 . Solid line: r = 2 pm; dashed 
line: r — 4 pm; dotted line: r — 6 pm; dot-dashed line: r = 8 pm; dashed double-dots line: r — 10 pm. 



r in the experiments of Dickinson et al. is unknown. Taking into account this uncertainty, the volume fraction of 
protein was systematically lowered while keeping the rest of the parameters fixed. On the one hand, it should be 
kept in mind that this procedure looks for an approximate form of the steric potential of casein and not for the exact 
value of the volume fraction. On the other hand, the calculation of the volume fraction requires knowledge of the 
intcrfacial area of the protein, its specific volume and the width of the protein layer. Hence, it summarizes the effect 
of several variables whose accurate values are unknown. Moreover, the introduction of a volume fraction as an input 
of the simulation avoids the dependence of the calculations on the actual value of the interfacial area of the /3-casein 
molecule (see the formula of V s t in Table I). Notice that the electrostatic potential used in the present simulations 
is also independent of the interfacial area of the protein, since it was parametrized using the charge density of the 
drops, not the effective charge of a protein molecule. As a result, a change in the area per molecule varies the charge 
of a surfactant molecule but does not change the value a, or the shape of the potential. 

A volume fraction of protein of <^> p = TO x 10~ 2 produces an inter-micellar potential similar in magnitude to the one 
reported by Tuinier and Kruif (Fig. 3). However, this volume fraction also produces large repulsive barriers between 
a drop and "the interface" (Fig. 5). 

Figure 6 show the potentials corresponding to the system composed by a micrometer size-drop and a 500 pm drop, 
when a volume fraction of protein equal to <p p = 4.7 x 10~ 3 is employed. The attractive van der Waals interaction 
prevails. Notice that the lowest potential curve corresponds to the biggest particle size. Despite this fact, the barrier 
height of the potential, AT^, increases with the particle size. Barriers of 2.4 kT, 3.8 kT, 6.1 kT, 8.2 kT and 11.0 kT are 
observed for drop sizes of r = 2 pm, r = 4 pm, r = 6 pm, r = 8 pm, r = 10 pm, respectively. 

The above potential was used in some preliminary evaluations of the coalescence time. These calculations are shown 
in the next section. As will be shown, the small barriers shown in Figure 6 still prevent the correct behavior of the 
coalescence time as a function of the particle radius. 
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FIG. 6: Total potential of interaction between a large droplet of 500 /jm, and a small droplet of micron-size. For this calculation 
the volume fraction of protein around the drops was assumed to be equal to <j> p — 4.7 x 10~ 3 . Solid line: r = 2/im; dashed 
line: r = 4 pm; dotted line: r = 6 /xm; dot-dashed line: r — 8 pm; dashed double-dots line: r — 10 fim. 



IV. EFFECT OF THE INTERACTION POTENTIAL, THE HYDRODYNAMIC FRICTION, AND THE 
BUOYANCY FORCE 

Table IV shows the results of preliminary calculations in which the potentials of Fig. 6 are employed. The average 
coalescence time between a drop and the interface was evaluated assuming non-deformable spherical droplets and the 
tensors of Stokes (outer region) and Taylor (inner zone). As shown in Fig. 7(a)-7(c), the coalescence time increases 
as a function of the particle radius due to the augment of the potential barrier with the particle size. 

Notice that in this case, the increase of the coalescence time shown in Figs. 7(a)-7(c) is independent of the initial 
distance of separation between the particle and the interface (d = 20 nm, 100 nm, 750 nm). It is also independent of 
the deformability of the drops, since a similar tendency is observed for the case of deformable droplets. This trend 
is contrary to the effect of the buoyancy force, which promotes a decrease of the coalescence time as a function of 
the particle radius. However it agrees with the fact that the diffusion tensor decreases with the size of the particles, 
increasing their coalescence time. Thus, the increase of the repulsive barrier reinforces the effect of the hydrodynamic 
friction and predominates over the buoyancy force. These results are in complete contradiction with the experimental 
evidence that shows a decrease of the coalescence time as a function of the particle radius for micrometric drops. 
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0.13751 1.52714 


7.10532 


14.69749 


14.38758 


Fig. 7(a) 


100 


0.32029 1.81169 


7.49218 


15.12774 


14.75408 


Fig. 7(b) 


750 


1.86863 2.71276 


8.09156 


15.75886 


15.15175 


Fig. 7(c) 



TABLE IV: Average coalescence time for the total potential of Fig. 6 




FIG. 7: Dependence of the average coalescence time as a function of the particle radius for an initial distance of: (a) d — 20 nm, 
(b) d = 100 nm, (c) d = 750 nm. The total potential corresponds to the one depicted in Fig. 6 (Table IV). 



13 



d (nm) 


t 2 (s) 
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Behaviour 


20 


0.00873 0.01636 0.02243 0.02758 0.03503 


Fig. 8(a) 


100 


0.20471 


0.26402 0.28301 0.27541 0.26570 


Fig. 8(b) 


750 


1.73811 


1.30137 1.02320 0.84397 0.72334 


Fig. 8(c) 



TABLE V: Average coalescence time for the total potential of Fig. 9 
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FIG. 8: Dependence of the average coalescence time as a function of the particle radius for an initial distance of: (a) d — 20 nm, 
(b) d — 100 nm, (c) d = 750 nm. The total potential corresponds to the one depicted in Fig. 9 (Table V). 



Table V and Figs. 8(a)-8(c) show similar calculations in which the potentials depicted in Fig. 9 are employed. These 
potentials were obtained, using a smaller volume fraction of protein (</> p = 7.9 x 1CP 4 ). This volume fraction produces 
a very small steric barrier (AV < 0.25 kT) between K-casein micelles. It also eliminates the repulsive barriers shown 
in Fig. 6 for the drop-interface system, generating a total potential that is more attractive as the particle radius 
increases. Such potentials favor shorter coalescence times for larger particle radii. Hence, in this case, the attractive 
potential reinforces the effect of the buoyancy force and opposes the hydrodynamic resistance. 

Curiously, it is found that the coalescence time (r) increases with the size of the drops for an initial distance of 
d = 20 nm,. This calculation was repeated twice: first eliminating the Brownian movement of the particles, and 
second eliminating all repulsive potentials between them. However, similar results were found in all three cases. This 
indicated that the variation of r with r 4 was due to the effect of the diffusion tensor. Table VI shows that the value of 
/corr for the Taylor tensor (Tensor 777 in Table II) decreases sensibly as a function of the particle radius, promoting 
longer coalescence times for larger particle radius. 



Tj (/im) Ahxn/rl 4h 2 rj/rl 4fe 3 r»/r^ 



2 0.0101 


0.0504 


0.3780 


4 0.0051 


0.0254 


0.1905 


6 0.0034 


0.0171 


0.1280 


8 0.0026 


0.0129 


0.0968 


10 0.0021 


0.0104 


0.0780 



TABLE VI: The ratio Ahm/rl for hi = 20 nm, h 2 = 100 nm and h 3 = 750 nm. 

As the initial distance of approach between the particles increases, the range of action of the attractive potential 
and the hydrodynamic friction increase. The diffusion tensor of Stokes depends linearly on the distance of approach, 
but inversely on the particle radius. According to Table VI the ratio D^/Dq increases in absolute magnitude as the 
initial distance of approach increases. The particles diffuse more rapidly through a longer distance, and differences 
in the magnitude of the attractive potential are experienced by the drops during a longer period of time. As a result 
of this phenomenon, an increase in the initial distance of approach changes the variation of the coalescence time as 
a function of the particle radius. The interaction potential and the buoyancy force dominate at d = 750 nm, and the 
initial trend exhibited at d = 20 nm is reversed (see Table V and Fig. 8(a)-8(b)). 

These results indicate that the dependence of the coalescence time as a function of the particle radius observed in 
the experiments requires: (a) the absence of repulsive barriers between the emerging drops and the interface, and (b) 
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a sufficiently long distance of approach [d > 750 nm). 

Figures 9 and 10 show the final potentials employed in the present simulations for spherical and deformable drops. 
The total potential corresponding to each particle radius is completely attractive in both cases. The corresponding 
repulsive barriers illustrated in Figs. 11 and 12 are counterbalanced by a large van der Waals attraction, see Fig. 9 
and 10. 
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FIG. 9: Total potential of interaction between a large spherical droplet of 500 /im and a small spherical droplet of micron-size. 
For this calculation the volume fraction of protein around the drops was assumed to be equal to (f> p = 7.9 x 10~ 4 . Solid line: 
r — 2/im; dashed line: r = 4 pm; dotted line: r = 6 /xm; dot-dashed line: r — 8 /jm; dashed double-dots line: r — 10 /nm. 










-5 




-10 










-15 








-20 


Z3 




> 


-25 




-30 




-35 




-40 



100 



200 300 
d (nm) 



500 



FIG. 10: Total potential of interaction between a large deformable droplet of 500 /im and a small deformable droplet of micron- 
size. For this calculation the volume fraction of protein around the drops was assumed to be equal to (j> p = 7.9 x 10~ 4 . Solid 
line: r = 2 /jm; dashed line: r — 4/im; dotted line: r = 6/im; dot-dashed line: r = 8 pm; dashed double-dots line: r = 10 /xm. 
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FIG. 11: Individual contributions to the total potential between a drop and a planar interface for a spherical droplet of 



r = 10 /im. The steric potential was calculated using a volume fraction of protein equal to 
electrostactic; dotted line: steric. 



= 7.9 x 10" 



Dashed line: 



Figure 13 shows the interaction potential between two drops of r = 10 /im. This potential was calculated using 
the same parameters that generate the potentials shown in Figs. 9 and 10 for the interaction between a drop and 



15 





25 




20 


p 


15 


*o 








"c 


10 






> 





5- 

- 

2 4 6 8 10 12 

d (nm) 

FIG. 12: Individual contributions to the total potential between a drop and a planar interface for a deformable droplet of 
r = 10 /im. The steric potential was calculated using a volume fraction of protein equal to (f> p = 7.9 x 1CP 4 . Dashed line: 
electrostactic; dotted line: steric; dot-dashed line: dilational; dashed double-dots line: bending. 



a planar interface. In the case of two deformable drops, the repulsive contributions produce a significant potential 
barrier, even at </> p = 7.9 x 10 -4 . Hence, care must be taken when extrapolating the experimental measurements of a 
drop/interface system to the coalescence behavior of two deformable drops. 
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FIG. 13: Total interaction potential between two micron-size drops with a radius of 10 /im. The parameters of the potentials 
are equal to the ones employed in the calculation of Fig. 9 and 10 for the drop/interface model. V Bt (cj> p = 7.91 x 10~ 4 ). Solid 
line: spherical drops; dashed line: deformable drops. 

Based on the results discussed above (Tables IV and V) , and the order of magnitude of the experimental coalescence 
times, several initial distances of approach were tested. For spherical particles d — 5 /im, 10 /im, 15 /im, 20 /im, and 
30 /im. In the case of deformable drops, the order of magnitude of the experimental coalescence times is achieved at 
considerably shorter distances of separation. In this case initial distances d = 50 nm, 65 nm and 100 nm were found to 
be convenient. The results presented in the following figures correspond to those distances of approach that showed 
the closest agreement with the experimental measurements. 

All calculations included the effect of the buoyancy force, Eq. (5). The time step of the simulations was kept large 
enough (At, = At x D a /r} = 1.0) in order to guarantee that At > mD /kT. For drops with l^tm < n < 10 /im, 
At* = 1 corresponds to values between 2.0 x 10" 6 s and 2.0 x 10~ 2 s, (mD /kT = 1.9 x 10 _7 s and 1.9 x 10~ 4 s, 
respectively). 




V. RESULTS AND DISCUSSION 



Figure 14 illustrates the results of the calculations for spherical drops. Notice that only six different particle 
radii were used in order to find the best combination of tensors. The separation distance which showed the closest 
agreement with the experimental data was d = 15 /im. This distance corresponds to the value of /ii n i in Eq. (6). 
Basheva et al. [7] suggested a simple procedure to estimate the distance h- m - u at which spherical droplets change their 
bulk velocity drastically due to the proximity of an interface. According to Eq. (6), a plot of the product rr^ vs. 
1/Vj should approach a straight line. Figure 15 shows the result of applying this procedure to the data of Dickinson 
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et al. regarding hexadecane drops covered by /3-casein. From the slope and the intercept of the curve, a value of 
/ii n i equal to d = 19.4 /jm was obtained. This value appears to be sufficiently close to 15/Ltm if it is considered that 
the procedure employed to implement the tensors in the simulations is not completely equivalent to the fitting of the 
analytical equation published in Ref. (4). 
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FIG. 14: Average coalescence time vs. droplet radius for spherical drops initially separated by a distance of d = 15 jj,m. Stars: 
experimental data [6]; double-dashed dot line: Tensor A; dashed double-dots line: Tensor B; dotted line: Tensor C; grey dashed 
line: Tensor D; solid line: Tensor E; dashed line: Tensor F; dot-dashed line: Tensor G; grey dot-dashed line: Tensor H . 

According to Figure 14, all the simulations which include combinations of Stokes and Taylor tensors (Tensors D, E, 
G and H) come close to the experimental data. In the scale of this figure, Tensors G and H overlap, as well as tensors 
D and E. It appears from these results that the mobility of the liquid near the interface is not a very significant factor 
at least for this specific case. A close look at the curves evidences that Tensor D shows the best behavior. Tensor F 
underestimates the coalescence time, while Tensor C overestimates it. The behavior of Tensor F is partially due to 
the use of the Taylor tensor outside its region of validity. In regard to Tensor C it should be noticed that the quality 
of the prediction of the tensors which include the correction of Honig et al. markedly depend on the value of the 
reference radii employed. The equation of Honig et al. was formerly deduced for two equal spheres. This explains 
why a reasonable agreement is obtained with Tensor C, while the predictions of Tensors A and B are totally mistaken. 
In the last two cases, the radius of reference was estimated as the average radius between the emerging drop and 
the fixed 500 ^m drop. As a result, i? lies far away from the radii of the coalescing drops, producing unreliable 
coalescence times. 
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FIG. 15: Fitting of Eq. (6) to the data of Dickinson et al. [6]. Stars: experimental data; solid line: Stokes- Taylor equation 

In Fig. 16 we show the behavior between average coalescence time and the droplet radius for the case of deformable 
drops. The value of h[ n i was in all simulations of deformable droplets was of the order of 10 nanometers, three orders 
of magnitude different from the one deduced by application of Eq. (6) to the experimental data (see Fig. 15). Notice 
that one of the curves shows a non-monotonic variation of the coalescence time (Tensor J). This is due to the correction 

(2) 

of the diffusion tensor of Reynolds when /cVr > 1-0. The program has the option to use the set of tensors available 
whatever their value, and also has the option to correct the calculated tensor using the expression of Honig et at, if 

(2) 

the absolute value of /cVr > 1.0. The fact that the correction was implemented in one of the simulations that used 
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the Reynolds tensor (Tensor J), indicates that analytical formula of this tensor might fall outside its range of validity 
when the scheme of Fig. 1 is implemented. 

The lowest monotonous curve in Fig. 16, includes the calculations ran with an initial distance of 50 nm. These 
comprise Tensors I, L, and M. No significant difference was observed between the simulations ran with the tensor of 
Danov et al. at e$ = 0.1 and 1.0. The highest monotonous curve include the calculations that employed an initial 
distance of d = 100 nm. The results from tensors I, K and M overlap in this case. 

In the case of Tensor K, all three mechanisms of coalescence previously described were tested. These included film 
drainage, and surface oscillations with two different forms of counting the time (ry). However, in the present case, 
the effect of the buoyancy force and the van der Waals attraction (deterministic forces) prevails. In the absence of a 
significant repulsive barrier, the drops coalesce once they enter the deformation zone. Hence, the evaluation of is 
irrelevant regardless of the method of counting the coalescing time. Notice also, that the use of large time-steps in the 
simulations might induce the artificial coalescence of the particles since they might produce values of Tjj considerably 
larger than ryrij- However, the use of a large arbitrary time (i.e. 10 s) as input of the simulations led to the same 
result. Coalescence occurred very rapidly as a consequence of deterministic forces. This was further confirmed by 
setting the Brownian contribution of Eq. (7) equal to zero. This is achieved using an auxiliary variable called Damper 
which multiplies the value of the random contributions. Notice that in the present simulations each particle moves 
individually with the equation of motion Eq. (7). Hence, the random deviates are ascribed to each coalescing particle 
regardless the formation of a doublet. However, it is probable that during the drainage of the film, a couple of 
dcformable drops might move as a unique entity. The program has the option to decrease the effect of the Brownian 
motion of each particle multiplying their random deviates by a real number lower than one. 

Tensors K and L were tested in order to consider the possible occurrence of deformation at higher distances of 
approach than the ones predicted by Eq. (12). The values of /i ini predicted by this equation do not surpass a few 
nanometers. However, as shown in Fig. 13 for the case of non-deformable particles, the presence of a planar interface 
promotes a substantial viscous friction at distances of a few micrometers. Hence, the program has the possibility to 
introduce an arbitrary value of h ln i as an input of the calculation, but in this case the formula for the evaluation of 
the film radius cannot be used, and its value also needs to be estimated beforehand. 

The methodology described for the simulation of deformable droplets is likely to be very useful in those cases in 
which the radius of the particles is larger, and a substantial repulsive potential exists. In the presence case, these 
procedures do not produce a better result than the simple combination of the tensors of Reynolds and Stokes (Tensor 
I). However, the fact that deformable droplets do not reproduce the experimental trend is not surprising. According 
to Ivanov et al. [14] the analytic expression of hi n i in terms of the interaction forces and the disjoining pressure does 
not have a solution for droplet radii smaller than 83 /im. This means that the surface of the drop remains convex and 
there is no formation of a film. The present result support that prediction for radii between 1 /im < < 10 /im. 
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FIG. 16: Average coalescence time vs. droplet radius for deformable drops. Stars: experimental data [6]; dashed line: Tensor 
J (d = 50 nm); dot-dashed line: Tensor K (d = 50 nm); dotted line: Tensor L (d = 50 nm); dashed double-dots line: Tensor I 
(d = 50 nm); double-dashed dot line: Tensor I (d = 100 nm). 

In order to make a closer comparison between the simulation and the experimental data, fifteen additional simula- 
tions corresponding to intermediate particle radii between 1 and 10 /im were run for the case of spherical particles. For 
these calculations, the two tensors that produced the closest agreement with the experimental points were employed 
(Tensor D and G in Fig. 14). Figure 17 shows the results of the calculations. The agreement between the simulations 
and the experimental data is very good for the case of tensor D. Notice that the standard deviation of the calculations 
is indicated in the figure. It is remarkable that the error bars increase monotonically with the decrease of the particle 
radii. Such dependence was observed by Dickinson using lysozyme as a surfactant, but the magnitude of the errors 
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was not plotted in Ref. [6] for the case of the /3-casein protein. It is likely that the error bars evidence the effect of the 
Brownian motion on the trajectory of the particles. However, a close analysis of the trajectories of each simulation 
was not carried out. 



FIG. 17: Best simulation results for the coalescence time of a micron-size drop pressed by buoyancy against a planar interface 
(d=15/im). Stars: experimental data [6]; solid line: Stokes- Taylor law (Tensor D); dashed line: Stokes- Taylor law (Tensor G). 
Error bars were approximated by the standard deviation of 1000 simulations 

The coalescence times predicted by Tensor G follow the curvature of the experimental points but do not get close 
enough. If the initial distance of separation is changed so that the experimental point corresponding to = 4 microns 
is reproduced, the rest of the coalescence times predicted lie above the experimental curve (for instance, r = 17.6 s 
for r, = 2 /mi, instead of 12.5 s as it was experimentally found). Moreover, the initial distance of separation has to 
be increased up to 30 fim in order to reproduce the point of r-j = 4 fim. This distance is considerably larger than 
the 19.2 ^m deduced by application of Eq. (6) to the experimental data of Dickinson ct al. [6]. Consequently, it 
appears that under the experimental conditions, the surface of the droplets acquire enough surfactant to behave as 
an immobile interface, despite the relatively short time of contact between the drops and the protein solution. 

VI. CONCLUSION 

In the absence of a significant repulsive barrier, the present simulations confirm the analytical predictions of Basheva 
et al. [7]. According to our calculations, the experimental behavior of hexadecane drops at a water/hexadecane 
interface [6] can be reproduced assuming spherical droplets that move with a combination of Stokes and Taylor 
tensors. The variation of the coalescence time as a function of the particle size, predicted by the simulations of 
deformable drops, does not follow the experimental trend. This suggests that micron size droplets behave as non 
deformable droplets in agreement with previous theoretical work [8] . 

In was also shown that the presence of a strong repulsive barrier between the emerging drop and the interface might 
completely change the behavior of the coalescence time as a function of the particle radius. In this regard, the initial 
distance of approach between the particle and the interface plays a significant role. 

Acknowledgements 

The authors acknowledge Dr. Aileen Loszan and Dr. Jhoan Toro-Mendoza for useful discussions and computational 
time. 



[1] K. P. Velikov, O. D. Velev, K. G. Marinova and G. N. Constantinides, J. Chem. Soc, Faraday Trans. 93, 2069 (1997). 

[2] R. R. Dagastine, R. Manica, S. L. Carnie, D. Y. C. Chan, G. W. Stevens and F. Grieser, Science 313, 210 (2006). 

[3] G. B. Webber, S. A. Edwards, G. W. Stevens, F. Grieser, R. R. Dagastine and D. Y. C. Chan, Soft Matter 4, 1270 (2008). 

[4] H. Casanova, J. Chen, E. Dickinson, B. S. Murray, P. V. Nelson and M. Whittle, Phys. Chem. Chem. Phys. 2, 3861 (2000). 

[5] S. R. Deshiik and K. D. Papadopoulos, J. Coll. Sci 174, 302 (1995). 

[6] E. Dickinson, B. S. Murray and G. Stainsby, J. Chem. Soc, Faraday Trans. 1 84, 871 (1988). 




50 



10 



r(um) 



19 



[7] E. S. Basheva, T. D. Gurkov, I. B. Ivanov, G. B. Bantchev and B. Campbell and R. P. Borwankar, Langmuir 15, 6754 
(1999). 

[8] I. B. Ivanov and P. A. Kralchevsky, Coll. Surf. A 128, 155 (1997). 

[9] G. I. Taylor, Proc. R. Soc. Lond. A 103, 58 (1923). 
[10] D. L. Ermak and J. A. McCammon, J. Chem. Phys. 69, 1352 (1978). 
[11] G. Urbina-Villalba, Int. J. Mol. Sci 10, 1 (2009). 

[12] J. Toro-Mendoza, A. Lozsan, M. Garcia-Sucre, A. J. Castellanos and G. Urbina-Villalba, Submitted to Physical Review E 
(2009). 

[13] T. G. M. van de Ven, Colloidal Hydrodynamics (Academic, Padstow, 1989), pp. 1-51. 
[14] I. B. Ivanov, K. D. Danov and P. A. Kralchevsky, Coll. Surf. A 152, 161 (1999). 

[15] K. D. Danov, N. D. Denkov, D. N. Petsev, I. B. Ivanov and R. Borwankar, Langmuir 9, 1731 (1993). 
[16] K. D. Danov, D. N. Petsev, N. D. Denkov and R. Borwankar, J. Chem. Phys. 99, 7179 (1993). 
[17] A. Scheludko, Adv. Coll. Int. Sci. 1, 391 (1967). 

[18] I. B. Ivanov, B. Radoev, E. Manev and A. Scheludko, Trans. Faraday Soc. 66, 1262 (1970). 
[19] E. D. Manev and A. V. Nguyen, Adv. Coll. Int. Sci. 114-115, 133 (2005). 
[20] E. D. Manev and J. K. Angarska, Coll. Surf. A 263, 250 (2005). 

[21] E. Dickinson and V. J. Pinfield and D. S. Home and F. A. M. Leermakers, J. Chem. Soc, Faraday Trans. 9, 1785 (1997). 
[22] A. V. Makienski. V. B. Fainerman, M. Bree, R. Wiistneck, J. Kragel and R. Miller, J. Phys. Chem. B 102, 417 (1998). 
[23] J. Maldonado-Valderrama, M. J. Galvez-Ruiz, A. Martin-Rodriguez and M.A. Cabrerizo-Vflchez, Langmuir 20, 6093 
(2004). 

[24] H. C. Hamaker, Physica IV, 1058 (1937). 

[25] A. Lozsan, M. Garcia-Sucre and G. Urbina-Villalba, Phys. Rev. E. 72, 061405 (2005). 

[26] A. Lozsan, M. Garcia-Sucre and G. Urbina-Villalba, J. Coll. Int. Sci 299, 366 (2006). 

[27] P. A. Kralchevsky, T. D. Gurkov and I. B. Ivanov, Coll. Surf. 56, 149 (1991). 

[28] P. A. Kralchevsky and T. D. Gurkov, Coll. Surf. 56, 101 (1991). 

[29] E. P. Honig, G. J. Roebersen and P. H. Wiersema, J. Coll. Int. Sci 36, 97 (1971). 

[30] P. Mills and P. Snabre, Europhys. Lett. 25, 651 (1994). 

[31] C. W. J. Beenakker and P. Mazur, Phys. Lett. A 91, 290 (1982). 

[32] C. W. J. Beenakker and P. Mazur, Physica A 126, 349 (1984). 

[33] G. Urbina-Villalba, M. Garcia-Sucre and J. Toro-Mendoza, Phys. Rev. E. 68, 061408 (2003). 

[34] T. D. Gurkov and E.S. Basheva, Encyclopedia of surface and colloid science, Vol. 4 (Taylor & Francis, 2002), p. 2773. 

[35] G. Urbina-Villalba, A. Lozsan, K. Rahn and M. S. Romero-Cano, Comput. Phys. Comm. 180, p. 2129 (2009). 

[36] A. Vrij, J. Coll. Int. Sci 19, 1 (1964). 

[37] A. Vrij, Discuss. Faraday Soc. 42, 23 (1966). 

[38] A. Vrij and J. Th. G. Overbeek, J. Am. Chem. Soc. 90, 3074 (1968). 
[39] D. Kashchiev and D. Exerowa, J. Coll. Int. Sci 77, 501 (1980). 
[40] A. Nikolova and D. Exerowa, Coll. Surf. A 149, 185 (1999). 

[41] S. Damodaran and C. S. Rao, Food Colloids: Fundamentals of Formulation (The Royal Society of Chemistry, 2001), p. 
165. 

[42] D. E. Graham and M. C. Phillips, J. Coll. Int. Sci 70, 403 (1979). 
[43] D. E. Graham and M. C. Phillips, J. Coll. Int. Sci 70, 415 (1979). 
[44] D. E. Graham and M. C. Phillips, J. Coll. Int. Sci 70, 427 (1979). 
[45] J. Chen, E. Dickinson and G. Ivesion, Food Structure 12, 135 (1993). 
[46] J. Chen and E. Dickinson, J. Sci. Foo Agric. 62, 283 (1993). 

[47] T. D. Dimitrova, F. Leal-Calderon, T. D. Gurkov and B. Campbell, Adv. Coll. Int. Sci. 108-109, 73 (2004). 
[48] E. Leclerc and P. Calmettes, Phys. Rev. Lett. 78, 150 (1997). 

[49] E. Dickinson and M. Golding and M. J. W. Povey, J. Coll. Int. Sci 185, 515 (1997). 
[50] R. Tuinier and C. G. de Kruif, J. Chem. Phys. 117, 1290 (2002). 
[51] C. G. de Kruif, J. Dairy Sci. 81, 3019 (1998). 

[52] H. -M. Wang , C. -Y. Li, S. -J. Chen, T. -W. Cheng and T. -L. Chen, J. Membr. Sci. 238, 1 (2004). 

[53] A. H. Nhaesi and A. -F. A. Asfour, J. Chem. Eng. Data 45, 991 (200). 

[54] D. Mobius and R. Miller, Proteins at Liquid Interfaces (Elsevier, 1998). 

[55] E. Dickinson, Int. Dairy J. 9, 305 (1999). 

[56] G. Urbina-Villalba, Langmuir 22, 5968 (2006). 

[57] J.N. Israelachvili, Intermolecular and surface forces (Academic Press, 1998), p. 90. 
[58] T. L. McMeekin, M. L. Groves and N. J. Hipp, J. Am. Chem. Soc 71, 3298 (1949). 
[59] E. Dickinson, D. S. Home, J. S. Phipps and R. M. Richardson, Langmuir 9, 242 (1993). 



20 



